In [99]:
Image.open("/home/xpwang/Downloads/colposcope3/train/benign/20130902133200.jpg")
Out[99]:
多阶段的分割算法¶
- Intensity $I = \frac{R + G + B}{3}$ Saturation $S = 1 - \frac{\min(R,G,B)}{I} $ $$ I > 0.4 * I_{max}\ \ and\ \ S < 0.6 * S_{max} $$
- gradient magnitude $ |\triangledown G| = \sqrt{(\frac{\partial I}{\partial x})^2 + (\frac{\partial I}{\partial y})^2}$ $$ |\triangledown G| > 0.15 $$
- Gaussin Mixture Model cluster
概率密度函数
$$
p(x_{i}) = \sum_{k=1}^{K}\pi_{k}P(x_i|j) \\
= \sum_{k=1}^{K}\pi_{k} \mathcal{N}(x_i|\mu_k,\Sigma_k)
$$
- S -- V feature space, where S and V are Saturation and Value of the standard HSV color space
- 似然函数： $$ \prod_{i=1}^{N}p(x_i) $$
- log似然函数： $$ \sum_{i=1}^{N}\log p(x_i) = \sum_{i=1}^{N}\log (\sum_{k=1}^{K}\pi_{k} \mathcal{N}(x_i|\mu_k,\Sigma_k)) $$
- EM算法：引入因变量z，似然函数变为：$$\prod_{i=1}^{N}p(x_i,z_i) = \prod_{i=1}^{N}p(z_i)p(x_i|z_i) = \prod_{i=1}^{N} \prod_{k=1}^{K}(\pi_k)^{z_k}\mathcal{N}(x_i|\mu_k,\Sigma_k)^{z_k}$$
- log似然：$$\sum_{i=1}^{N}\sum_{i=1}^{N}z_{ik}\{\log \pi_k + \log \mathcal{N}(x_i|\mu_k,\Sigma_k)\}$$
- 估计$z_k$,根据后验概率取期望：$$z_k = 0 * p(z_k = 0|x) + 1 * p(z_k = 1 | x) \\ = p(z_k = 1 | x) \\= \frac{p(x)p(x|z_k =1)}{p(z_k)} \\= \frac{\pi_k\mathcal{N}(x|\mu_k,\Sigma_k)}{\sum_{j=1}^{K}\pi_j\mathcal{N}(x|\mu_j,\Sigma_j)}$$
mixture of Gaussian the EM is the following.
Repeat until converged:
E-step. Compute the expectations, or responsibilities $r_{j,i}$, of the latent variables $v_i$, where $j$ represents the classes, $y$ the data points, and $i$ the $ith$ points: $$\begin{equation} r_{j,i} = \frac{\omega_j P(y_i~|~\mu_j,\Sigma_j)}{\sum_j \omega_j P(y_i~|~\mu_j,\Sigma_j)} \end{equation} $$
M-step. Compute the maximum likelihood parameters given these responsibilities:
(1) The weigths for the mixing proportions, $\omega_j$, where $N$ is the number of data points: The estimate of the probability of a class as an outcome is the proportion of times it appears: $$\begin{equation} \omega_j ← \arg\max \sum_i \sum_j r_{j,i} log P(z_i = j) = \sum_i \sum_j r_{j,i} log \omega_j = \frac{1}{N}\sum_i r_{j,i} \end{equation} $$
(2) The means and the variances for the Gauss distributions, $\mu_j$, $\sigma_j$ : We want to maximize this term: $$\begin{equation} \arg\max \sum_i r_{j,i} log P(y_i~|~v_i = j) = \sum_i r_{j,i} log Gaussian(y_i; \mu_j, \sigma_j) \end{equation} $$ $$\begin{equation} \mu_j ← \frac{1}{\sum_i r_{j,i}} \sum_i r_{j,i} y_i \end{equation} $$ $$\begin{equation} \sigma_j ← \frac{1}{\sum_i r_{j,i}} \sum_i r_{j,i} (y_i - \mu_j)^2 \end{equation} $$
填充算法¶
In [6]:
def multile_stage(image,I_lamda = 0.7,S_lamda = 0.1) :
# 转换成array
arr = np.array(image)
# 计算intensity 和 sataration
arr_I = np.mean(arr,axis =2)
arr_S = 1 - np.min(arr,axis=2) / (arr_I+1e-7)
# 根据阈值生成mask
ise = arr_I > I_lamda * np.max(arr_I)
sse = arr_S < S_lamda * np.max(arr_S)
# 计算梯度
magnitude = generic_gradient_magnitude(arr_I,sobel)
mse = magnitude > np.percentile(magnitude,97)
#第一阶段的效果
select_1 = ise & sse
print(ise.sum(),sse.sum(),mse.sum())
arr_1 = arr.copy()
arr_1[:,:,0][select_1] = 0
arr_1[:,:,1][select_1] = 0
arr_1[:,:,2][select_1] = 0
# 第二阶段的效果
select = ise & sse | mse
arr[:,:,0][select] = 0
arr[:,:,1][select] = 0
arr[:,:,2][select] = 0
return arr_1,arr,select_1,select,mse
In [7]:
def gassin_process(image,mask):
arr = np.array(image)
mask = mask.copy()
arr_I = np.mean(arr,axis =2)
arr_S = 1 - np.min(arr,axis=2) / (arr_I+1e-7)
hsv = color.rgb2hsv(arr)
S = arr_S[mask]
V = hsv[:, :, 2][mask]
X = np.vstack((S, V)).T
model = GaussianMixture(n_components=2, covariance_type="full", max_iter=20000)
model.fit(X)
resoult = model.predict(X)
mask[mask][resoult == 0] = False
arr[:, :, 0][mask] = 0.0
arr[:, :, 1][mask] = 0.0
arr[:, :, 2][mask] = 0.0
return arr,mask
In [8]:
def fill_process(image,mask,kenel_size = 20) :
image = np.array(image)
mask = mask.copy()
mask = dilation(mask,selem=square(5))
image[:, :, 0][mask != 0] = 0.0
image[:, :, 1][mask != 0] = 0.0
image[:, :, 2][mask != 0] = 0.0
while(mask.sum() != 0) :
init = mask.sum()
shape = mask.shape
print(mask.sum())
for i in range(mask.shape[0]) :
for j in range(mask.shape[1]) :
nz = 0
sum_nz = np.array([0,0,0])
if mask[i,j] != 0 :
for ni in np.array(range(kenel_size)) - kenel_size/2:
x = int(i + ni)
for nj in np.array(range(kenel_size))- kenel_size/2 :
y = int(j + nj)
if x >0 and x < shape[0] and y > 0 and y < shape[1] and image[x,y,:].sum() != 0 and image[x,y,:].sum() < 245*3:
nz += 1
sum_nz += image[x,y,:]
if nz != 0 :
image[i,j,:] = sum_nz / nz
mask[i,j] = False
# print(mask.sum())
if init == mask.sum() :
# image[:, :, 0][self.select] = np.median(self.arr[:,:,0])
# self.arr[:, :, 1][self.select] = np.median(self.arr[:,:,1])
# self.arr[:, :, 2][self.select] = np.median(self.arr[:,:,2])
break
return image
In [10]:
plt.figure(figsize=(40,40))
plt.subplot(1,2,1)
plt.imshow(image)
plt.title("Image")
plt.subplot(1,2,2)
plt.imshow(test1)
plt.title("after Intensity and Saturation selection")
Out[10]:
In [12]:
plt.figure(figsize=(40,40))
plt.subplot(1,2,1)
plt.imshow(test2)
plt.title("after gradient magtitude selection")
plt.subplot(1,2,2)
plt.imshow(test3)
plt.title("after GMM")
Out[12]:
In [13]:
plt.figure(figsize=(10,10))
# plt.subplot(1,2,1)
plt.imshow(test4)
plt.title("after fill process")
Out[13]:
第二种分割算法¶
- stage1:
- $c_E = 0.2989 · c_R + 0.5870·c_G + 0.1140·c_B$
- $r_{GE} = \frac{P_{95}(c_{G})}{P_{95}(c_E)}$ $r_{BE} = \frac{P_{95}(c_{B})}{P_{95}(c_E)}$
- $c_{G}(x_0) > r_{GE} * T_1 \mid c_{B}(x_0) > r_{BE} * T_1 \mid c_E(x_0) > T_1 $
- stage2:
- $\varepsilon(x_0) = \max \{ \frac{c_R(x_0)}{c^*_R(x_0)},\frac{c_G(x_0)}{c^*_G(x_0)},\frac{c_B(x_0)}{c^*_B(x_0)} \} $
- $c^*(x) 是一个经过填充，滤波后的图片$
- $\varepsilon(x) > T_2$
- stage3:
- $\frac{1}{N}\sum_{i=1}{N}|grad(E_n)| > T_3 \land N > N_{min}$
In [95]:
contour_list = find_contours(modual_2(),0.1)
plt.figure()
fig,ax = plt.subplots(figsize=(10,10))
ax.imshow(modual_2(),cmap = "gray")
for contour in contour_list :
ax.plot(contour[:,1],contour[:,0])
In [94]:
image
Out[94]:
比较好的方法¶
- 算法要有很好的鲁棒性，不能在一张图片上过拟合
- 思想就是既要考虑召回率又要考虑准确度，加一点，减一点
- 算法要考虑到数据本身的实际情况，不能简单的拿来套用，比如要考虑消除图片上的tag，图片反光的特征
In [40]:
def sr_mark(image) :
arr = np.array(image)
sele_var = np.var(arr,axis=2)
# sele_min = np.min(arr,axis = 2)
choose = sele_var.copy()
choose[:,:] = False
choose[50:120,80:570][sele_var[50:120,80:570] <100] = True
choose[120:,80:643][sele_var[120:,80:643] <100] = True
# choose[60:633,63:575][sele_min[60:633,63:575] > 180] = True
arr[:,:,0][choose != 0] = 0
arr[:,:,1][choose != 0] = 0
arr[:,:,2][choose != 0] = 0
# arr[:,:,0] = median_filter(arr[:,:,0],size=5)
# arr[:,:,1] = median_filter(arr[:,:,1],size=5)
# arr[:,:,2] = median_filter(arr[:,:,2],size=5)
choose = dilation(choose,selem=square(5))
return arr,choose
In [41]:
def fill_process(image,mask,kenel_size = 20) :
mask = mask.copy()
image = np.array(image)
image[:, :, 0][mask != 0] = 0.0
image[:, :, 1][mask !=0] = 0.0
image[:, :, 2][mask != 0] = 0.0
while(mask.sum() != 0) :
init = mask.sum()
shape = mask.shape
for i in range(mask.shape[0]) :
for j in range(mask.shape[1]) :
nz = 0
sum_nz = np.array([0,0,0])
if mask[i,j] != 0 :
for ni in np.array(range(kenel_size)) - kenel_size/2:
x = int(i + ni)
for nj in np.array(range(kenel_size))- kenel_size/2 :
y = int(j + nj)
if x >0 and x < shape[0] and y > 0 and y < shape[1] and image[x,y,:].sum() != 0 and image[x,y,:].sum() < 245*3:
nz += 1
sum_nz += image[x,y,:]
if nz != 0 :
image[i,j,:] = sum_nz / nz
mask[i,j] = False
# print(mask.sum())
if init == mask.sum() :
# image[:, :, 0][self.select] = np.median(self.arr[:,:,0])
# self.arr[:, :, 1][self.select] = np.median(self.arr[:,:,1])
# self.arr[:, :, 2][self.select] = np.median(self.arr[:,:,2])
break
return image
In [48]:
def background(image,back,percent = 99.5) :
back = back.copy()
image = image.copy()
back[:,:,0] = median_filter(back[:,:,0],size=10)
back[:,:,1] = median_filter(back[:,:,1],size =10)
back[:,:,2] = median_filter(back[:,:,2],size= 10)
mask = np.min(image / back ,axis =2 )
image[mask > np.percentile(mask,percent),:] = back[mask > np.percentile(mask,percent),:]
return image,mask,back
In [ ]:
_,mask = sr_mark(image)
paint_2 = fill_process(image,mask)
paint_4,mask_4,back_4 = background(paint_2,paint_2)
In [15]:
plt.figure(figsize=(8,8))
plt.imshow(paint_2)
plt.title("after var remove and fill")
Out[15]:
In [215]:
plt.figure(figsize=(8,8))
plt.imshow(paint_4)
plt.title("after background fillprocess")
Out[215]:
In [82]:
show(1)
In [57]:
show(2)
In [58]:
show(3)
In [79]:
show(4)
In [80]:
show(5)
In [62]:
show(9)
In [64]:
show(np.random.randint(1)%80)
In [72]:
show(np.random.randint(80))
In [78]:
show(np.random.randint(80))
参考文章¶
- Arnold M, Ghosh A, Ameling S, et al. Automatic segmentation and inpainting of specular highlights for endoscopic imaging[J]. EURASIP Journal on Image and Video Processing, 2010, 2010(1): 814319.
- Zimmerman-Moreno G, Greenspan H. Automatic detection of specular reflections in uterine cervix images[J]. Medical Imaging, 2006, 6144(1): 61446E-61446E.